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Abstract: Reconstruction algorithms for diffuse optical tomography (DOT) 
and bioluminescence tomography (BLT) have been developed based on 
diffusion theory. The algorithms numerically solve the diffusion equation 
using the finite element method. The direct measurements of the 
uncalibrated light fluence rates by a camera are used for the reconstructions. 
The DOT is self-calibrated by using all possible pairs of transmission 
images obtained with external sources along with the relative values of the 
simulated data and the calculated Jacobian. The reconstruction is done in the 
relative domain with the cancelation of any geometrical or optical factors. 
The transmission measurements for the DOT are used for calibrating the 
bioluminescence measurements at each wavelength and then a normalized 
system of equations is built up which is self-calibrated for the BLT. The 
algorithms have been applied to a three dimensional model of the mouse 
(MOBY) segmented into tissue regions which are assumed to have uniform 
optical properties. The DOT uses the direct method for calculating the 
Jacobian. The BLT uses a reduced space of eigenvectors of the Green's 
function with iterative shrinking of the permissible source region. The 
reconstruction results of the DOT and BLT algorithms show good 
agreement with the actual values when using either absolute or relative data. 
Even a small calibration error causes significant degradation of the 
reconstructions based on absolute data. 
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1. Introduction 

Bioluminescence imaging (BLI) is a noninvasive technique that can monitor biological 
processes at a molecular level and visualize disease progression and response to treatment [1- 
4]. In its simplest form, a sensitive camera is used to image light that reaches the surface from 
the internal bioluminescence source. These images can provide useful information but they 
are not quantitative. Bioluminescence tomography (BLT) attempts to reconstruct the three 
dimensional distribution of the bioluminescent source (i.e. optical power per unit volume) and 
to provide true quantitative information about its magnitude and location [5-9]. Many 
different strategies have been proposed for BLT [10-18] but, in general, three components are 
necessary: a forward model of light propagation from the internal source to the animal 
surface, a calibration of the camera that links the image (e.g. counts per pixel) to the emittance 
(optical power per unit area), and a reconstruction algorithm, usually iterative, that estimates 
the source distribution that minimizes the difference between the observed images and those 
calculated by application of the forward model. 

The second component, camera calibration, typically requires three sub-components: an 
assumption that the emission from the animal surface is Lambertian [19], modeling of the 
camera optics that transfers the light to the image plane, and measurement of the absolute 
camera response on a pixel-by-pixel basis. Recent work by Kienle and Foschum showed that 
the Lambertian approximation is a poor description of the true surface emission [20]. In 
addition, measurement of the camera response and modeling of the optics are challenging 
technical problems. In this paper we propose an alternative strategy for BLT. In addition to 
the bioluminescence image, a transmission image (or images) is acquired using an external 
fiber optic source (or sources) of known optical power. The reconstruction algorithm uses the 
ratio of the internal source image to the external source image as input data. Because the 
angular dependence of the emission at any point on the surface is very likely the same for both 
sources, and because the same camera is used to capture both images, the camera calibration 
problem is reduced to a simple mapping of the animal surface to the camera image plane. In 
our system [21] the surface anatomy is acquired by an integral x-ray CT scanner but other 
methods, such as structured light imaging, have been applied [22]. 

Several investigators have concluded that better BLT reconstructions can be achieved if 
forward calculations take into account the true three dimensional distribution of scattering and 
absorption coefficients [23,24]. In previous papers [25,26] we have described a strategy 
whereby this information can be obtained by diffuse optical tomography (DOT). Transmission 
images acquired with external fiber optics sources and the BLI camera can be used to 
reconstruct the optical properties which, in turn, are used in the BLT algorithm. Our previous 
algorithms assumed that a perfect camera calibration (as described in the previous paragraph) 
was available to link the transmission images to the actual emittance produced by each source. 
In this paper we propose an alternative strategy for DOT: instead of "absolute data", ratio 
images or "relative data" obtained with different source positions are used. If, for example, 
nine external sources are employed and nine transmission images are acquired, 36 ratio 
images can be formed and used as input to a DOT algorithm. The relative power of each of 
the fiber optic sources must also be known. In the simulations presented in this paper we 
demonstrate that the optical property reconstructions based on these ratios are generally not as 
accurate as those based on absolute calibration. However, the BLT solutions based on the 
relative data are comparable in accuracy to those based on absolute data. In addition, we show 
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by example that the relative method outperforms the absolute method if the calibration 
process is imperfect. 

2. DOT algorithm 

Although we will use only the continuous wave (CW) solution, for generality we begin with 
the frequency domain diffusion equation and the Robin-type boundary condition at a specific 
wavelength given by [27,28] 

-V-*(M)V + f^(M) + i®j *(r,a>,X) = s(r,a>,A.) (retl), (1) 

[l + 2 K -(M)^(r)// a (M)]^(r,A,«) = 0 (re SO), (2) 

where (p{r,X,co) is the light fluence rate at position r, wavelength 1, and light modulation 
frequency / given by a> = 2nf ; O and 5Q are the domain and its boundary, respectively, and 
the source distribution is given by s(r,co,X). The spatial distribution of the tissue optical 
properties at wavelength A is given by the absorption coefficient /u a (r,/L) and the diffusion 
coefficient k(t,X), where the diffusion coefficient is defined by s(r, X) = 

[3(// a (r,/l) + /^(r,/l))] , where fj.' s (r,X) is the reduced scattering coefficient 

fi[ = jU s (l-g), Ms is tne scattering coefficient, and g is the anisotropy factor. «(r)is a unit 
vector pointed outwardly normal to <3Q, and f is derived from Fresnel's law as [29] 
£ = ((2 / (1 - R 0 )) -1 + |cos 6 C | 3 ) / (l - |cos 6> c | 2 ) , where the critical angle G c = sin" 1 (\ln) and 

R 0 = (n-lf / (n + lf , n is the tissue refractive index, and the refractive index of the 

surrounding medium is assumed to be 1 . 

Equation (1) can be solved numerically using the method of finite elements (FE) [30-32], 
and the discrete version of Eq. (1) is given by 

A{X)*{X) = s{X), (3) 

where A is the forward matrix obtained by discretizing the diffusion operator in Eq. (1). 
Assuming the object can be segmented into N regions where each element within the region 
has the same scattering and absorption coefficients, the number of unknown optical 
coefficients is 2N and the optical properties vector can be written in its transpose form as 
jU T =[ju' sl ,...,ju' sN ,ju al ,...,ju aN ], where /u' m and n an are the scattering and absorption 
coefficients of region n. Differentiating both sides of Eq. (3) with respect to /u n gives 

dA d(j> ds ... 
<j> + A—!- = = 0, (4) 

where the differentiation of the source with respect to the optical properties gives zero, since 
the source is not a function of the optical properties. The derivative of the forward matrix A 
with respect to //„ is obtained by changing /u n by ±A// n and calculating the difference in the 
forward matrix AA . The first order derivative of the light fluence rate with respect to the 
optical properties is obtained from Eq. (4): 

y- = -A-*™-*. (5) 
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The first order derivative of the light fluence rate with respect to the optical properties at 
the detectors around the object is given by 



(6) 



where d is the detector index and </> d is a vector of length N d , where N d is the number of 

detectors and the matrix A' 1 (of,:) consists of the d rows and all columns of the matrix A' 1 . If 

the detectors do not coincide with the FE nodes at the boundary, interpolation can be used to 
get the value from the surrounding nodes. It is assumed that the values of the light fluence rate 
at the detector locations can be obtained from the measurements of the corresponding pixels 

of the CCD camera. For a phasor representation of the light fluence rate where <j) d = (p d e' e " , 
Eq. (6) can be written as 



i<p d e' 



3"„ 



3"„ 



dfi n dn„ dn n dn n 



(7) 



Therefore 



31og(%) 



f 



= real 



dju n 



= imag 



(8) 



The light fluence rate data captured by the detectors around the object can be 
approximated by a first order Taylor expansion such that 



log( % ) = log( %0 ) + X^^«, 



0„ 



8 + V^L?„ 



(9) 



where (p d0 and 0 do are the magnitude and phase of the light fluence rate corresponding to the 
initial guess jU 0 . In matrix form, Eq. (9) is given by 



log(<2> rf )Wlog(<z> rf0 ) 



Slog(^) 31og(^ rf ) 



3//, 



dju 2 



(10) 



which can be concisely written as 



-Jd/j. 



(11) 



For multiple external sources, measured data are obtained for each source and so Eq. (11) 
for N. sources is modified to 
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V 7 ^ J 



djj, 



(12) 



or, concisely, <S> d = <£> d0 +Jd/u. This equation can be solved iteratively for d jU such that 



d M =(j T jy'(j T (® d -® d0 )). 



(13) 



Equation (13) can be used to reconstruct the optical properties for absolute (or completely 
calibrated) data, where the magnitude and phase of the light fluence rate at the object 
boundary can be obtained by heterodyne detection. In the case of non-contact camera 
detection, usually only the magnitude is measured (i.e.,/= 0). 

If the measured light fluence rate is not calibrated, the magnitude of the light fluence rate 
at the detectors is multiplied by an unknown vector of calibration coefficients a and the phase 
is shifted by an unknown vector of phase delay ft. Therefore the measured light fluence rate 
<j> m is given by 



a<p d e 







^og(<p d f 


+ 








V °d J 


v P ) 



log (or) 



(14) 



and for multiple source measurements 



3)> 



log(«) V 
P 



r log(«) 

p 



(15) 



A system of equations to reconstruct the optical properties as in Eq. (12) can be obtained 
by using the relative uncalibrated measured data to remove the magnification factor vector a 
and the phase delay vector /? such that 



<x>! 



-d) 1 

m 



O • - (J) 1 

m n 

<J> 3 -<J) 2 

m m 

-<J} 2 



-<t 2 



V m m J 
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<t> 3 -d) 2 
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f -Wi ) 




J 3 — J \ 








J 3 — J 2 
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J 4 -J 2 




J If, ~ J 2 
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dju, (16) 
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or, concisely, <S m = <S rf0 + Jd/u. 

Therefore, the optical properties can be reconstructed iteratively from the relative 
measured data O m and the relative simulated data <D rf0 at the initial guess ju 0 by using the 
following equation: 

d M = (j T jf(j T (^ m -& d0 )). (17) 

Equation (17) can be used to reconstruct the optical properties by using the relative 
measurement, and the reconstructed optical properties at iteration i are given by 

fi t =//,_, +dfi v (18) 

The relative simulated data <D rf0 are recalculated at every iteration, using the updated 

optical properties //. . The relative Jacobian J can be updated at every iteration if there is 

high nonlinearity or less frequently if there are no large changes in the optical properties 
during these iterations. The solution is obtained when the maximum relative error between the 

relative measured data O m and the relative simulated data O rf0 goes below a certain limit. 
3. BLT algorithm 

The diffusion equation in Eq. (3) can be written as [26] 

cp(X) = G(X)s(Z), (19) 

where G(/l) is the Green's function matrix at wavelength X and s(X) is the bioluminescence 

source distribution. Assuming that the source magnitude is not a function of wavelength or, 
equivalently, the source spectrum is known in the range of measurement, the light fluence rate 
at the detectors on the object surface is given by 

(p(d;X) = G(d,R;X)s(R) or <p d (X) = G dR (X)s R , (20) 

where d is the detector index at the boundary <3Q , and R is the permissible source region 
index which contains the index of nodes of the finite element mesh inside the domain Q.. 

The light fluence rate measured by the CCD camera is related to the light fluence rate at 
the "detectors" on DO. by 

<P m =a-9 d => <P m =a-G dR s R , (21) 

where a is an unknown scaling vector that represents all calibrations required to obtain the 
right values of the light fluence rate at the detectors. Using Eq. (21), the logarithm of the 
measured light fluence rates corresponding to the bioluminescence source s R and the external 
sources s lt s 2 ,..., s N are given by 

log p s m * = log a + log ( G dR s R ) 
log <P'2 =loga + log(G rf s,) 

log qf^ = log a + log(G d s 2 ) (22) 
log^f' =\oga + \og( < G d s N ) 
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where G d =G(d, ) is the Green's function where all columns of the matrix are included, 

which means the permissible region is the whole object. The external sources s u s 2 ,..., s N 

used for the DOT reconstruction and the measurements resulted from these sources are also 
used for calibrating the bioluminescence measurements. The calibrating vector a can be 
removed from Eq. (21), and a calibrated equation can be obtained as 



G dR s R = exp 



logfl? " — Z( lQ g^ -M G A)) I => G dR s R = ft , (23) 



where <p"* are the calibrated bioluminescence measurements. For multi-spectral 
measurements, Eq. (23) for N wavelengths in a normalized form can be given by 



G d M)' 



max 



G«(^ 2 )/max 



^(A 2 )/max 



9 s : (A iV )/max(^(A JV )) 



(24) 



or, concisely, G dR s R = q> m . A symmetric system of equations can be obtained from Eq. (24): 

{G T dR G dR )s R =G T dR cp m . (25) 

As described in detail in a previous paper [33], Eq. (25) can be solved by assuming that 
the bioluminescence source s R is written as a superposition of the eigenvectors of the matrix 

[GrftfGjs] corresponding to non-zero eigenvalues such that s R ='^ j "_ x a i v i , where 

v = [vj,v 2 v M ] is a set of eigenvectors of the matrix [G^G^I corresponding to non-zero 

eigenvalues and M is the number of chosen eigenvectors. The vector 
a = \v T (G^G^vJ [v^G^)^,! consists of the coefficients of the expansion [33]. The 

equation described above is solved several times for different permissible source regions 
starting from the whole object size and ending with a few nodes. For every solution, the 

objective function f=\G dR [:,R)s R -<p m \ , which gives the difference between the calculated 

and measured light fluence rate, is calculated. After all iterations, the solution that minimizes 
the objective function is considered the best solution. The shrinking of the permissible region 
after every iteration is done by removing source nodes from the domain R corresponding to 
low source magnitude as described in [33] 

4. Results and discussion 

To demonstrate these algorithms, a highly detailed anatomy for a laboratory mouse (MOBY) 
[34] is used to create the FE mesh. The 3D object (Fig. 1(a)) represents a finite element mesh 
of a section of a mouse abdomen where adipose, liver and spleen, bowel and intestine, kidney, 
and bone tissues are identified and assigned realistic wavelength-dependent optical properties. 
It is assumed that all elements within each region have the same optical properties. There are 
9 optical fibers used as external sources located on the ventral surface of the mouse as shown 
in Fig. 1(b). It is assumed that these libers produce Gaussian sources located one transport 
length inside the object (~1 mm) with a maximum value of 1 nW/mm 3 and standard deviation 
of 1 mm. The measurements corresponding to these sources are used for the DOT 
reconstruction and BLT calibration at all wavelengths. The detectors shown in Fig. 1(b) are 
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2183 triangular surface elements and form a field of view of 225 degrees on the dorsal 
surface. It is assumed that the light fluence rate at these detectors can be obtained from CCD 
camera images by appropriate geometrical mapping [35]. 



Spleen 




(a) (b) 

Fig. 1 . (a) A 3D finite element mesh for the object used in the simulation. The mesh is a regular 
Delaunay tetrahedral with 17966 nodes and 99140 tetrahedrons generated by iso2mesh [36]. 
The object represents a section of a mouse abdomen which illustrates 5 different tissue regions: 
(1) adipose, (2) liver and spleen, (3) bowel and intestine, (4) kidneys, and (5) bone, (b) A 
schematic of the object showing the sources and detectors set up and the coordinate system. 

In the first step, the optical properties of all tissue regions at the wavelengths 580, 590, 
600, 610, 620, and 630 nm were reconstructed using the DOT algorithm described above. The 
light fluence rate at each detector is determined for excitation of each source fiber. However 
these measurements are not calibrated and are shifted in both magnitude and phase from the 
correct values <D d as shown in Eq. (15). Equations (16) and (17) show that the optical 
properties can be reconstructed by using relative measured data, relative simulated data, and 
the relative Jacobian as shown in Eq. (16). The use of relative values in Eq. (16) cancels out 
the shift in magnitude and phase of the measured data from the correct calibrated data shown 
in Eq. (15). The DOT algorithm described above can be used for both frequency domain 
/ > 0 or continuous wave (CW) / = 0 . In case of CW, the phase 9 and the derivatives 

^^/qjj m Eq- (10) will be zeros, which can be removed from the equations or left, as there is 

no effect on the results. The initial estimates used for scattering and absorption for all tissues 
at all wavelengths are 1 mm" and 0.01 mm" respectively. The DOT algorithm stops when 
the maximum relative error of the measured and simulated light fluence rates is less than 1% 
or the number of iterations exceeds 60. 

Figure 2 shows a comparison of the relative percentage error in reconstructing the optical 
properties for CW at wavelength 590 nm using absolute data, relative data, and imperfect 
calibrated data. The forward model, Eq. (3), is used to simulate the measured data using the 
actual optical properties and 2% Gaussian noise is added to the calculated data to simulate 
realistic conditions. The true scattering and absorption coefficients were calculated using the 
data reported in [24,37] and are shown in Tables 1 and 2 respectively. The corresponding 
relative errors of reconstructing scattering, absorption, and effective attenuation 

ju e/f = ^3 ( ju' s + ju a )/u a are shown in the figure. The errors using relative data are larger than 
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Fig. 2. The relative percent error in the reconstructions of scattering fi , absorption // ,and 
effective attenuation ju , for the 5 tissue regions (1) adipose, (2) liver, (3) bowel, (4) kidney, 

and (5) bone at 590 nm using absolute, relative, and imperfect calibrations. Note the split 
vertical scale. 

when using absolute data. The maximum relative errors obtained in the range 580 nm to 630 
nm are 2% for scattering, 3% for absorption, and 1% for effective attenuation when using 
absolute data while the maximum errors are 9% for scattering, 12% for absorption, and 4% for 
effective attenuation when using relative data. The simulation took 32 minutes for 60 
iterations on a pc with Intel processor with a clock of 3.33 GHz and 24 GB RAM. The 
Jacobian was recalculated using the new reconstructed optical properties every 10 iterations. 

In Fig. 2 we also show an example of the results obtained when the calibration process is 
imperfect. There are many possible sources of error in this multi-step process but here we use 
a simple illustration. One way to calibrate the camera is to use a reference object with known 
optical properties [38]. We have assumed that this reference object has the same external 
shape as the Moby phantom and uniform internal optical properties. We further assume that 
the reduced scattering coefficient is nominally 1.00 mm -1 but actually 1.05 mm -1 and the 
absorption coefficient is nominally 0.0100 mm -1 but actually 0.0105 mm -1 . We performed 
forward calculations of the transmission for each source-detector pair under nominal and 
actual conditions. The ratio is the imperfect calibration factor used in all subsequent 
simulations including the BLT calculations. This 5% error in scattering and absorption 
coefficients is relatively small and likely underestimates the cumulative error in a real 
calibration process. With this imperfect calibration the maximum relative errors are 8% for 
scattering, 55% for absorption, and 35% for effective attenuation. 

For BLT, to simulate the measured data q>'* using Eq. (23), the true source distribution s R 
and the Green's function G dR calculated using the true optical properties at each wavelength 
given in Tables 1 and 2 were used. 2% Gaussian noise was added to simulate realistic data: 
(p s * in Eq. (23). Two cases are simulated for <p^ in Eq. (23): the first case is to set each 

Table 1. True [24,37] reduced scattering coefficients (mm 4 ) for each region in the 
heterogeneous object at six different wavelengths 





580 nm 


590 nm 


600 nm 


610 nm 


620 nm 


630 nm 


Adipose 


1.30 


1.29 


1.28 


1.27 


1.26 


1.25 


Liver 


0.79 


0.77 


0.76 


0.75 


0.74 


0.72 


Bowel 


1.37 


1.35 


1.32 


1.29 


1.27 


1.24 


Kidney 


2.80 


2.73 


2.66 


2.60 


2.53 


2.47 


Bone 


3.08 


3.01 


2.93 


2.86 


2.80 


2.73 
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Table 2. True [24,37] absorption coefficients (mm ') for each region in the heterogeneous 
object at six different wavelengths 





580 nm 


590 nm 


600 nm 


610 nm 


620 nm 


630 nm 


Adipose 


0.0073 


0.0043 


0.0021 


0.0014 


0.0010 


0.0008 


Liver 


0.6474 


0.3993 


0.1899 


0.1200 


0.0824 


0.0647 


Bowel 


0.0198 


0.0128 


0.0063 


0.0040 


0.0028 


0.0023 


Kidney 


0.1209 


0.0746 


0.0356 


0.0226 


0.0156 


0.0123 


Bone 


0.1040 


0.0670 


0.0325 


0.0207 


0.0142 


0.0112 



element in a equal to 1 in Eq. (22) when the data are calibrated and match the correct values 
at the detectors around the object. In this case there is no need to use external sources for 
calibrations and ^ = qf* in Eq. (23). The second case is to set each element in a equal to 
1000 in Eq. (22) to simulate uncalibrated data and normalization with external sources is 
needed as described in Eq. (23). To simulate the measured data due to external sources (p s m , 
the Green's functions G d calculated using the true optical properties at each wavelength were 
used in Eq. (22). The external sources are simulated by Gaussian sources as described before 
and the elements of a equal 1000. In Eq. (23), the Green's functions G dR and G d were 

calculated using the reconstructed optical properties, by the DOT algorithm described above, 
using absolute data, relative data, and imperfect calibrated data. The permissible region R was 
initially chosen to be a cylinder of 20 mm diameter and height 20 mm with its center located 
at (0, 3, 3). At each iteration, the permissible region is slightly shrunk by an arbitrary factor y 
The permissible region initially containing 5010 mesh nodes is finally shrunk to 10 nodes in 
60 iterations, y = (5010/10) <1/60> = 1.1092. The 60 iterations took approximately 30 minutes on 
a PC with Intel processor with a clock of 3.33 GHz and 24 GB RAM. The execution time 
varies greatly according to the choice of mesh size, the volume of the initial permissible 
region, and number of iterations. It is useful to start with a small size mesh and fewer 
iterations to get a quick initial solution before solving the complete problem. This can provide 
a more focused permissible region and hence greatly reduce the simulation time and give 
more accurate results. 

Figure 3 shows transverse and coronal cross sections for the actual and reconstructed 
source distributions. The transverse section passes through the two kidneys at z = 3 mm and 
the coronal section passes through the two kidneys at y = 3 mm . The actual sources are 

spheres of 4 mm diameters located in the two kidneys at (-5, 3, 3) and (6, 3, 3) with 
magnitudes of 2 nW/ mm 3 and 1 nW/ mm 3 respectively and total spatially integrated powers 
of 67.2 nW and 39 nW respectively. 

Two objective metrics have been calculated to assess the accuracy of the BLT 
reconstruction of each source. The first is simply the total power integrated over three 

dimensions, ^s(i)V (z) , where V(i) is the volume of the element i in the FE mesh and the 

summation includes all elements that include the source. The second we define as the 
"normalized magnitude error" calculated according to 

Z actual s -\ rcconst / • \ \ j /- / -\ i V ' actual / -xts/ - \ 
\s (t)-s (i)\V(i) / 2_s (j)V(j)- 

i i 

This element-by-element summation is sensitive to the shape and position of the 
reconstructed source and is a more stringent metric than the total power. 

The reconstructed sources in Figs. 3(c) and 3(d) using absolute data and in 3(e) and 3(f) 
using relative data show good agreement with the actual sources in Figs. 3(a) and 3(b) in 
terms of location, total power reconstructed (within 0.75%) and magnitudes with normalized 
magnitude errors of 0.27. There are no significant differences in the reconstruction results 
when using absolute or relative data. The reconstructed sources using imperfect calibrated 



#172525 - $15.00 USD Received 13 Jul 2012; revised 3 Oct 2012; accepted 9 Oct 2012; published 1 1 Oct 2012 
(C) 2012 OSA 1 November 2012 / Vol. 3, No. 1 1 / BIOMEDICAL OPTICS EXPRESS 2804 




Source 1 Source 2 
Total power 67.2 nW 39.0 nW 

Magnitude 2 nW/mm' I nW/mnr 




Total power 67.4 nW 39.0 nW 

Normalized q.24 0 . 04 
magnitude error 



(c) 



(d) 




Total power 67.7 nW 39.1 nW 



Normalized 
magnitude error 



(1.2(1 



(e) 



(f) 




Total power 7 1 .2 n W 40.5 nW 



Normalized 



magnitude error 




Fig. 3. BLT reconstruction of two uniform bioluminescence sources localized in the left and 
right kidneys (4 mm diameters). The "anatomy" is shown in Fig. 1. The actual sources in 
transverse section (a) and coronal section (b) have magnitudes 2 nW/ mm 3 and 1 nW/ mm 3 for 
the left and right source respectively. The reconstructed sources are shown in (c) and (d) using 
absolute data, (e) and (f) using relative data, and (g) and (h) using imperfect calibrated data. 
The total powers of the actual and reconstructed sources and the normalized magnitude errors 
of the reconstructed sources using different reconstruction methods are shown in the figure. 
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Source I Source 2 
Total power 308.3 nW 316.3 nW 

Magnitude 1 nWW I nW/mnr 



Total power 309.5 nW 315.1 nW 





Total power 328.6 nW 332.4 nW 

Normalized () , , () , 5 

magnitude error 



0 0.5 1 1.5 o 0 5 1 15 

(9) (h) 

Fig. 4. BLT reconstruction of large distributed sources that fill the entire two kidneys. The 
magnitude of the actual sources in (a) and (b) is 1 nW/mm , and the reconstructed sources are 
shown in (c) and (d) using absolute data, (e) and (f) using relative data, and (g) and (h) using 
imperfect calibrated data. The total powers of the actual and reconstructed sources and the 
normalized magnitude errors of the reconstructed sources using different reconstruction 
methods are shown in the figure. 
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(g) (h) 

Fig. 5. BLT reconstruction of a large non-uniform source distribution that fills the entire right 
kidney. The actual source in (a) and (b) has a hot spot of 4 mm diameter with magnitude 2 
nW/mm 3 which is twice the magnitude of the background. The reconstructed source in (c) and 
(d) using the absolute data, (e) and (f) using the relative data, and (g) and (g) using the 
imperfect calibrated data. The total powers of the actual and reconstructed sources and the 
normalized magnitude errors of the reconstructed sources using different reconstruction 
methods are shown in the figure. 
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data are shown in Figs. 3(g) and 3(h) with total power reconstructed (within 6%) and 
normalized magnitude error of 0.54. 

Figure 4 shows the result of the reconstruction of large uniform sources that fill the entire 
two kidneys; the actual sources in Figs. 4(a) and 4(b) are of uniform magnitude 1 nW/mm 3 . 
The total powers of the sources are 308.3 nW in the left kidney and 316.3 nW in the right 
kidney. The reconstructed sources in Figs. 4(c) and 4(d) using absolute data and in 4(e) and 
4(f) using relative data show comparable results in terms of the total power reconstructed 
which is within 0.63% of the actual values and the normalized magnitude errors which are 
0.14 and 0.17 respectively. The reconstructed sources using imperfect calibration in Figs. 4(g) 
and 4(h) have total power reconstructed within 6.3% and normalized magnitude error of 0.54. 

Figure 5 shows the BLT reconstruction for a large non-uniform source that fills the entire 
right kidney. There is a 4 mm diameter "hot spot" of 2 nW/mm 3 located at (6, 3, 3) inside the 
kidney while the rest of the organ "background" of 1 nW/mm 3 as shown in Figs. 5(a) and 
5(b). The total power of the actual source is 355.3 nW which is reconstructed within 1% when 
using absolute data in Figs. 5(c) and 5(d) or relative data in 5(e) and 5(f). The hot spot can be 
recognized in both cases but the location is slightly more accurate when absolute data are 
used. The normalized magnitude error is 0.2 using absolute data and 0.25 while using relative 
data. The total power is reconstructed within 5.3% using imperfect calibration and the 
normalized magnitude error is 0.39. The hot spot location is less accurate than that using 
absolute or relative data. 

5. Conclusion 

Self-calibrated algorithms for three dimensional DOT and BLT based on the finite element 
solution of the diffusion equation have been developed. Both algorithms avoid absolute 
calibration of the detectors as by using ratios, camera images can be used for reconstructing 
the optical properties and bioluminescence sources. The algorithms have been applied to 
simulations based on a highly detailed anatomy for a laboratory mouse (MOBY) and show 
good agreement between the actual values of optical properties and bioluminescence sources 
and the reconstructions. The DOT algorithm used the direct method to calculate the Jacobian. 
Using relative values for the measured data, simulated data, and Jacobian showed good 
reconstruction results compared to using absolute values. The same measurements due to 
external sources used for DOT at each wavelength were also used to calibrate the 
bioluminescence measurements. The reconstruction results for BLT using relative calibration 
show good agreements with the results using absolute data. The BLT algorithms used the 
eigenvectors expansion of the Green's function as a reduced space for the reconstruction, and 
iterative shrinking of the permissible source region. Small multiple sources and large 
homogeneous and inhomogeneous source distributions were reconstructed using both absolute 
and relative data. The results showed good agreement between the reconstructed and actual 
sources in terms of location and total power and magnitude. The BLT reconstructions based 
on relative data are slightly inferior to those based on absolute data obtained with perfect 
calibration. However, any calibration procedure will be imperfect, and we have demonstrated 
that even small errors in the calibration can lead to performance that is worse than that based 
on relative data. 
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